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Abstract 



Two previous papers in this series have presented a study of the growth of 

52 ' hadronic bubbles during the cosmological Quark-Hadron transition, treating 

the material within each phase as a single perfect fluid. Here, we extend the 

^ ' analysis to include the effects of long-range energy and momentum transfer 

c^ \ by weakly and electromagnetically interacting particles. After a short review 

of the formalism adopted, we discuss the numerical strategies used in the 

computer code which has been constructed in order to solve this system 

of equations. Results for the growth of single hadronic bubbles are also 

presented. 
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I. Introduction 

At about 10 /US after the big bang it is thought that strongly interacting 
matter underwent a transition from a plasma of free quarks and gluons to one in 
which the quarks were confined within hadrons (predominantly pions). Current 
lattice-gauge calculations [1] favour the view that this change may well have been 
a continuous one, but the possibility of it being a first order phase transition is 
by no means ruled out and considerable interest continues to be focussed on the 
astrophysical consequences which might have arisen if it was, indeed, first order 
(see, for example [2-4]). 

The work described here (which is part of an ongoing research programme [5- 
8]) is within the scenario of a first-order transition starting with the nucleation of 
hadronic bubbles in a slightly supercooled quark-gluon plasma. The bubbles then 
proceed to grow (with the quark-gluon plasma being progressively transformed into 
a hadronic one) until eventually they coalesce and give rise to disconnected quark 
regions which then shrink and probably disappear completely (possibly leaving 
behind a significant baryon-number inhomogeneity) [9] but might instead reach a 
stable configuration composed of strange quark matter [10]. 

In the earlier papers of the series, attention has been focussed on studying the 
hydrodynamics of the growth of single hadronic bubbles during the initial stages 
of the transition where it makes sense to consider the material in each phase as a 
perfect fluid composed only of the strongly interacting matter. While the transition 
involves only these particles in a direct way, an important role is also played by 
other particles present, which can interact with them through the electromagnetic 
and weak interactions: primarily photons, electrons, muons and their antiparticles 
(electromagnetic and weak interactions) and neutrinos and antineutrinos (weak 
interaction only). All of these have mean free paths long compared with that 
of the strongly interacting matter and can provide a mechanism for long-range 
transport of energy and momentum through the strongly interacting fluid. Since 
the leptons concerned are essentially massless, both they and the photons can be 
treated as components of a generalized "radiation fluid" and the problem is then 
one of relativistic radiative transfer. 

During bubble growth, the effect of this transport becomes signiflcant when 
the radius Rg of the bubble surface becomes roughly comparable with the mean 
free path A of the particles concerned (~ 10"^ fermi for the electromagnetic in- 
teraction and ~ 1 cm for the weak interaction). When Rg <^ A, the bubble is 



essentially transparent to the radiation, which can then be ignored, while when 
Rs ^ X the coupling is essentially complete on relevant length-scales so that the 
radiation and strongly interacting matter move together as a single fluid. Clearly, 
the process of coupling can, in principle, occur twice during the bubble growth, 
but because the behaviour is similar in each case, we will discuss here only the 
one occurring between the strongly interacting and electromagnetically interacting 
particles. Identical considerations apply also for the coupling with the neutrinos, 
the only difference, apart from the scale, being the different number of degrees of 
freedom into which the energy liberated by the transition is channeled. 

This paper is concerned with adding the effects of radiative transfer to our 
earlier scheme of calculation for following bubble growth [7], thus extending its 
range of validity. Subsequently, we aim to apply a similar scheme to the last stage 
of the transition (shrinking away of the disconnected quark regions) which is likely 
to give rise to the most interesting consequences from an astrophysical point of 
view. The formalism used for handling the radiative transfer has been the subject 
of an earlier paper [8] and only a brief outline of this will be presented again here. 

The plan of the paper is as follows. In Section II we give a short review 
of the equations and then Section III describes the outline of our computational 
scheme which uses a combined characteristic and Lagrangian finite-difference ap- 
proach. Some particular problems arise in the computational implementation and 
we explain how these have been overcome by means of introducing new variables 
for the radiation fiuid. Section IV concerns the setting up of initial conditions. 
Section V describes tests carried out and results obtained and Section VI is the 
conclusion. Throughout, we use a system of units in which c = h = kg = 1 and a 
space-like signature (— , +, +, +). Greek indices are taken to run from to 3 and 
Latin indices from 1 to 3; covariant derivatives are denoted with a semi-colon and 
partial derivatives with a comma. 

II. Relativistic Hydrodynamical Equations 

This section contains a short review of the equations used in the rest of the 
paper. A full discussion of these has already been presented in [8] to which the 
reader is referred for further details and references. 

The transfer of energy and momentum is considered as taking place between 
a standard fluid, composed of the strongly interacting particles and any other 



particles effectively moving together with them, and a radiation fluid consisting of 
those particles responsible for the long-range transfer of energy and momentum. 
We treat the transfer using the PSTF (Projected Symmetric Trace Free) tensor 
formalism of Thorne [11], [12] in which the relativistic radiative transfer equation 
is rewritten in terms of an infinite hierarchy of differential equations involving 
PSTF moments which is then truncated by introducing suitable closure relations. 
For spherical symmetry, the tensor formalism becomes effectively a scalar one. We 
made the truncation at second order and obtained 
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(1) 

(2) 

W2 = /b^o, (3) 

where wq is the energy density of the radiation, t^i is the radiation flux and W2 is 
the anisotropy scalar of the radiation (all measured in the local rest frame of the 
standard fluid); sq and si are energy and momentum source functions and a and 
b are metric coefficients of the spherically symmetric line element 

ds^ = -a^dt^ + h'^d^i'^ + R^l.de'^ + sin^O rfv?^), (4) 

where /U is a comoving radial coordinate and R is the associated Eulerean coordi- 
nate (the Schwarzschild circumference coordinate). 

The quantity fj^ appearing in the closure relation (3) is a variable Eddington 
factor which can take values ranging from 0, for complete isotropy (which could 
be caused by the medium being extremely optically thick), to 2/3 for complete 
anisotropy (which might arise when the medium is very optically thin). The 
expression used for it has to be arrived at on the basis of physical considerations 
and for our case we have used 
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where A is the effective mean-free-path of the radiation particles. In our picture, 
the bubble is first nucleated within a uniform and isotropic quark medium and the 
radiation field (which interacts with the quark medium on suitably large scales 
and is in thermal equilibrium with it) initially shares these properties {i.e. wq is 



constant everywhere and wi and W2 are zero). The radius of the newly-formed 
bubble is small compared with A and the medium is essentially transparent to the 
radiation on this scale. As the bubble starts to expand, the radiation quantities 
deviate from their values at the time of nucleation primarily as a result of the 
Doppler effect arising from the motion of the standard fluid rest frames with re- 
spect to that of the radiation field (which is, so far, remaining uniform in its own 
frame). These Doppler corrections can be calculated analytically and the results 
were presented in the Appendix of [7] . Solely on the basis of this consideration, 
one finds that /^ is given by (8/9)w^(l + 4w^/3)~^. As the bubble grows to dimen- 
sions comparable with A, there is progressive coupling between the radiation and 
the standard fiuid on the relevant length-scales and this interaction tends in the 
direction of making the radiation more isotropic as seen from the standard fiuid. 
This effect is approximated by multiplying the Doppler term in (5) by a correction 
factor (in the large brackets) which has the effect of producing the right behaviour 
in the optically thin and optically thick limits and giving a physically plausible 
join between them. 

The appropriate value of A is not known accurately but, on the basis of 
elementary considerations, we have taken A ~ 10"^ fm. Tests made in order to 
investigate the sensitivity of the numerical code to the values adopted for the 
various parameters will be discussed in Section V. 

For the source moments sq and si we use the expressions 

so = j{e-wo) + {so)^^, (6) 

.1 = -^, (7) 

where (so)sc i^ ^ term expressing the contribution due to scatterings and e is 
the energy density for radiation in thermal equilibrium with the standard fiuid. 
Assuming that it roughly follows a black-body law, e can be written as 

. = 9,g)T^ (8) 

with g^ being the number of degrees of freedom of the radiation fiuid and T^ the 
temperature of the standard fluid. 

Obtaining a suitable expression for (so)sc ^^ ^^^^ straightforward. While 
detailed derivations have been made for simpler applications [13-14], the lack of 



precise knowledge about the interaction processes in the present case, has led us 
express (so)sc by the simple absorption and emission term 

(«o)sc = yle-w'o), (9) 

where 0.2 is an adjustable coefficient ranging between zero and one. Fortunately, 
the results of the numerical calculations turn out not to depend sensitively on the 
value chosen; a discussion of this will be given in Section V. 

Equations (l)-(3), describing the processes of transfer between the standard 
fluid and the radiation fluid, need to be solved together with the hydrodynamic 
equations 
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M,^ = 4nR'Rje + wo + -wil (14) 

1 + "-^ =^fl„. (15) 

where p is the relative compression factor (which plays the same role as played 
by the rest-mass density in a non-relativistic fluid), F is the general relativistic 
analogue of the Lorentz factor, and w is the specific enthalpy {w = (e +p)/p)- 
The generalized mass function M can also be calculated using the alternative 
equation 

M^t = -'iTrR^R^t(p+i:Wo + -wi + W2y (17) 

Equations of state are required for both phases of the strongly interacting matter. 
The hadronic medium is taken to consist of massless, pointlike pions for which 
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while the quark phase is described by the M.I.T. Bag Model equation of state [15] 
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where B is the "Bag" constant. We take g^ = 37, (7/^ = 3 and these values need to 
be incremented by the relevant number of additional degrees of freedom g^ when 
any of the radiative particles are completely coupled to the strongly interacting 
matter. 

For treating the phase interface, we again use a characteristic scheme (as in 
[7]). The characteristic form of (1), (2), (10) and (12) is 
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which are to be solved along the forward and backward characteristics of the 
standard fluid dp = ± {a/b)cs dt (here Cg = (dp/de)^''^ is the local sound speed in 
the standard fluid), and 
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(21) 



which are to be solved along the forward and backward characteristics of the 
radiation fluid d/j = ± (a/6) (1/3 + fj^Y^'^dt , and where for compactness we have 
defined 

K = U n^^"""/,, • (22) 

u[l + Aw^ /6)Wi 

Supplementary equations calculating p, R and M are solved along the flow- 
lines of the standard fluid (i.e. advective directions d^ = 0) and are 

dp-^^dp~ — rft = 0, (23) 
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dR = audt, (24) 
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The configuration of characteristic curves adjacent to the interface is shown 
in Figure 1 for evolution of the system from time level t to a subsequent time level 
t + At. The dashed lines represent the forward and backward characteristics for 
the radiation fluid r, the full narrow lines are the equivalent characteristics for the 
standard fluids f, the vertical dotted line is a flow-line of the standard fluid in the 
quark phase and the heavy line is the worldline of the interface. 



Figure 1. The configuration of characteristic curves near the phase interface drawn 
in the Lagrangian coordinate frame. 

Note that the difference between the characteristic directions results from 
the difference between the sound speeds in the radiation fluid ((1/3 + IeY^"^) ^-i^d 
in the standard fluid (cg). If the latter were not relativistic, this difference would 
be large but in the present case Cg -^ 1/a/3 and the difference between the sound 
speeds is frequently very small. This leads to some considerable complications in 
the numerical solution of the equations which we will discuss in Section IV. 

To complete the solution at the interface, it is also necessary to introduce 
junction conditions across it and those for the energy and momentum of the stan- 
dard fluid can be obtained using the Gauss-Codazzi formalism [16, 17]. Taking 



the surface tension a to be independent of temperature, these junction conditions 
are 



[{e + p)ab]^ = 0, 
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(26) 
(27) 
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where [A]"^ = A+ — A~ , {A}"^ = A+ + A~, //^ is the interface location, /i^ = 
dfig/dt, f = (a^ — 6^/i^)^/^ and the superscripts ^ indicate quantities immediately 
ahead of and behind the interface [5] . 

Up to the time of the complete coupling, it is reasonable to neglect any 
interaction of the radiation fluid with the matter in the phase interface and so the 
energy and momentum junction conditions for the radiation are just continuity 
conditions: 
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Other supplementary junction conditions follow from continuity across the 
interface of the metric quantities R, dR/dt, and ds 
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and from the time evolution of the mass function M 
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At the time of bubble nucleation, conditions are essentially Newtonian so 
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[M]± = 471 Rla. 



(34) 
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One further equation is needed in order to complete the solution at the 
interface and this is an expression for the rate at which the quark matter passes 
across it. A suitable expression is obtained by setting the hydrodynamical flux F^ 
into the hadron region equal to the net thermal flux F^ into it 

where cti is an accommodation coefficient (0 < cti < 1). 

This completes the set of equations. In the next Section we will present the 
details of the computational scheme and discuss the techniques and the strategies 
used. 



III. Numerical Computations of Bubble Growth 

In order to solve the equations of the previous section, we have constructed a 
computer code following the lines of the one developed previously for calculations 
of the initial stages of bubble growth (see [7] for a full description) . As before, this 
employs a composite numerical strategy in which a standard Lagrangian finite- 
difference method is used to solve the hydrodynamical equations in the bulk of 
each phase, while the system of characteristic equations and junction conditions 
is solved for the grid zones immediately ahead of and behind the phase interface 
(which is tracked continuously through the grid). The radiation quantities wq and 
W2 are taken as "mid-zone" quantities while the radiation fiux wi is taken as a 
zone boundary quantity and calculated at the half time level (as for the velocity 
u - see [7] for details). 

The equations presented in Section II are general in nature and can be applied 
to a variety of situations. Normally, there would be no problem in doing this but 
some particular difficulties have arisen when applying them to the present case of 
bubble growth at the cosmological Q-H phase transition. Here, direct use of the 
radiation equations in the form given above leads to rapidly-growing instabilities 
which destroy the solution. After a series of experiments it was found that the 
difficulty originates in the very small deviation of wq and t^i away from their 
initial values during the early part of the bubble growth, and in the fact that 
the characteristic sound speed in the radiation fluid (1/3 + fj^Y becomes very 
close to the sound speed in the standard fluid Cg when the radiation is nearly 
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isotropic in the rest frame of the standard fluid {i.e. when /^ ^0). These 
features lead to production of cancellation errors in the solution of Eqs. (1) and 
(2) and near divergences in the characteristic form of the equations (21) (the 
expression (c^ — 1/3 — /^j ) appears in the denominator of several terms). Note that 
the near equality of the sound speeds only arises when the standard fluid is, itself, 
relativistic (with Cs ~ l/\/3). Also, it is a peculiarity of the present situation that, 
initially, the radiation is nearly isotropic in the rest frame of the standard fluid 
not because the medium is optically thick on the scale of the bubble but, rather, 
because of the assumed isotropy of the universe. 

For overcoming the cancellation errors, we have introduced new radiation 
variables deflned as the difference between the energy density, the flux and the 
shear of the radiation fluid and some reference values (indicated below by the 
superscript *) with the aim of performing an analytic cancellation of large terms 
in the equations leaving behind smaller "difference" terms. It turns out to be 
convenient to take these reference values to be those which would be measured if 
the only effect were that resulting from the motion of the fluid relative to a uniform 
radiation fleld having an energy density equal to that at the time of nucleation 
of the bubble, (wo)iv (These are the pure Doppler values mentioned earlier and 
calculated in the Appendix of [7].) Using a tilde to denote the new variables, we 
have 



wo = wo- {wo)* =wo- ( 1 + -w^ j (tyo)jv, (36) 

4 
wi =wi- (wi)* = wi + -ur{wo)^, (37) 

Q 

W2=W2-{w2)*=W2--u'^{wo)^, (38) 

y 



and equations (1) and (2) can then be rewritten as 
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where the partial derivatives with respect to ^ have been replaced by the equivalent 
derivatives with respect to R {i.e. d/dR = {Att R'^ p /V)d / d y) . Equations (39) and 
(40) are the new radiation hydro dynamical equations for the bulk of each phase; 
once the "tilde" variables have been computed, the values of wo,wi,W2 can be 
calculated from (36)-(38). Note that in (39), (40) the radiation variables which 
are multiplied by G are not transformed according to (36)-(38). This has been 
done to keep the expressions in a simpler form and because the contribution of 
these terms is small under the present circumstances. 

Using the new variables, the radiation characteristic equations become 
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The new form of the radiation junction conditions (28) and (29) is 
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Note that the characteristic equation no longer has terms with (c^ — 1/3 — J^,) 
in the denominator but it does have terms containing the ratio Wi/Z^ and these 
still give rise to numerical instabilities. However, this can be countered by further 
rewriting the equations in a form in which Wi/Z^ only appears as the coefficient 
of expressions which are small when /^ ^ 0. The central point in our strategy 
consists in isolating the group of terms which appears on the left-hand-side of the 
fluid characteristic equations (20) (and hence tends to zero when ij; — > and the 
fluid and radiation characteristics coincide). Details of the manipulation involved 
are given in the Appendix. This group of terms can then be conveniently handled 
using the differences between parameter values at the feet of the fluid and radiation 
characteristics. The revised form of the radiation characteristic equations (which 
is the one actually implemented in the code) is 
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where BT is the same as in (44). Note that the last two Unes of (47) are the terms 
on the left-hand-side of the standard fluid characteristic equations (20) . 



IV. Initial conditions 

Making use of the experience gained in the earlier work [7], the setting of 
initial conditions has been quite straightforward. We started with a single super- 
cooled hadronic bubble nucleated in mechanical and thermal equilibrium with its 
surroundings at a temperature T^ slightly below the critical temperature for the 
transition T^ . The equilibrium is an unstable one and any perturbation (continu- 
ing expansion of the universe, for example) will cause it to start growing. However, 
this growth is extremely slow and, in practice, it is not easy to follow with our 
code as numerical noise rapidly dominates. Our strategy then was to introduce 
a small artificial perturbation, decreasing the fiuid temperature inside the bubble 
by a small amount below its equilibrium value, and analytically tracing the effect 
of this on related quantities (including the velocity field of the standard fiuid). 
This successfully produced a suitable data set for starting the time evolution in a 
smooth and consistent way (see [7]). 

Since, at this stage, the radius of the bubble is small compared with A, the 
radiation fiuid is not significantly affected by the thermal perturbation and re- 
mains uniform and isotropic in its own frame. The initial conditions for wo,wi, W2 
(measured in the rest frames of the standard fiuid) are then those calculated from 
the Doppler formulae discussed earlier. 



14 



V. Tests and results 

As usual in numerical computations, the construction of the computer code 
was followed by a series of tests to eliminate errors and verify that the strategies 
used were satisfactory. One important test consisted in turning off the source 
functions and checking that the computed values of the radiation variables agreed 
with the analytical Doppler expressions. This revealed the problems discussed 
earlier. When these had been satisfactorily solved, the source functions were then 
turned on again and complete runs of the code were carried out. As the radius 
of the bubble increased (leading to increased coupling between the radiation and 
the standard fluid on relevant length scales) care was required as increasingly 
steep gradients of t^o appeared in the vicinity of the interface prior to complete 
coupling. Since structure on a scale smaller than the grid spacing can obviously 
not be resolved, it is necessary to be ready to switch on complete coupling in the 
equations at the appropriate moment. Some experimentation was required in order 
to do this in the best way. When this had been done, further tests were carried 
out in order to examine the sensitivity of the results to changes in the physical 
parameters and assumptions. In the following, we will first present results for a 
set of "canonical" parameter values and then discuss the effect of varying some of 
these. 

In Figs. 2 - 5, we show results from a run with T^ = 150 MeV, T^^/T^ = 0.98, 
(Jo = 1, tti = a2 = 1 and A = 10^ fm. Here ao = o'/T^ is a parameter commonly 
used to measure the relative strength of the surface tension. The value which we 
are taking for this is larger than currently preferred ones but we give results for this 
case to allow direct comparison with those of [7] . Figs. 2 and 3 show the behaviour 
of the velocity and energy density of the standard fiuid, at various times during 
the bubble growth, while Figs. 4 and 5 show the corresponding behaviour of the 
radiation energy density and fiux. All of these figures should be viewed together 
and it is useful, also, to make comparison with Figs. 2 and 4 of [7] which are for 
the equivalent calculation without inclusion of the radiation particles. (Note that 
for convenience in drawing these figures, the values of the variables at the centre 
of the bubble have been plotted at logioR{fm) = rather than at R = 0.) 

During the first part of the bubble expansion (Rg < 10^ fm), the standard 
fiuid variables behave in an identical way to that seen previously in the calcula- 
tion with no radiative transfer: the velocity of the interface progressively increases 
and a compression wave is pushed out into the surrounding quark medium. The 
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velocity profile in the quark phase is approximately solenoidal {u oc 1/-R^). The 
radiation variables at this stage have profiles which are almost exactly the Doppler 
ones produced by the motion of the fluid relative to an essentially uniform radia- 
tion field. In the previous calculations, where radiative transfer was not included, 
the standard fluid variables tended towards a similarity solution which was effec- 
tively attained for Rg > 10"^ fm. In the present calculation, the coupling together 
of the radiation and the standard fluid (which can be clearly seen in Fig. 4) starts 
to be effective before the former similarity solution is fully reached, causing flrst a 
distortion of the velocity proflle and subsequently a decrease in the peak velocity. 
(If a smaller value is used for ctq, both bubble nucleation and the attainment of the 
similarity solution occur at smaller values of Rs-) As the coupling becomes more 
complete on the relevant length scales (when R^ ~ IC^ fm), the peak of the radi- 
ation flux proflle becomes very narrow (the main part of the flux is concentrated 
exactly at the interface) and when it is no longer possible to resolve this on the 
grid we switch to total coupling. This involves setting to zero the radiation flux 
wi, the Eddington factor /^ and the source functions sq and si and augmenting 
the number of degrees of freedom for the standard fluid at the interface to include 
also those of the coupled radiation. The behaviour of the radiation at the inter- 
face is then included together with that of the standard fluid. Following the total 
coupling, the variables tend rapidly to a new similarity solution characterized by 
smaller velocities (the front now has to push a medium having larger inertia) and 
a smaller temperature jump across the interface. 



Figures 2-3. Velocity of the standard fluid u and energy density e. Different curves 
are for different times during the bubble growth; the dashed lines represent the 
values at the initial time. Here cto = 1 and A = 10 fm. 



Figures 4-5. Radiation fluid energy density Wq and energy flux w-^ . Here a2 = 1 
and A = 10'' fm. 
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Next, we turn to a discussion of the effect on the results of varying the values 
taken for some of the physical parameters. Bearing in mind the discussion already 
given in [7], we will concentrate here just on the mean-free-path for the radiation 
particles A (whose value has an uncertainty of about an order of magnitude) , and 
on the coefficient 02 appearing in the expression for the source moment sq (Eq. 
(9)). 



Figures 6-7. Radiation energy density w„ and standard fluid velocity u just ahead 
of the phase interface and just behind it (higher and lower curves respectively). The 
curves are the result of calculations with different values of A; R^ is the bubble 
radius and Ofg = 1. The curves are plotted so that the hadronic region is on the left 
of the discontinuity tracing the phase interface. 

The effect produced by varying A is illustrated in Figures 6 and 7, which show 
the behaviour during the bubble expansion of the values of wq and u measured just 
ahead of the interface and just behind it. The different values of A used (which 
are represented with the different styles of line) are 5 x 10^ fni, 10"^ fni and 10^ fm, 
proceeding from left to right. As one would expect, use of a larger value of A has 
the effect of causing the coupling to occur when the bubble has reached a larger 
value of i?o. 



Figures 8-9. Radiation energy density Wq and standard fluid velocity u just ahead of 
the phase interface and just behind it. The curves are the result of calculations with 
different values of the non-conservative scattering coefficient aj. Here A = 10'' fm. 

Figures 8 and 9 show the effect of varying the coefficient for the non- 
conservative scatterings (02)- The different values of 0^2 used are 0, 0.5 and 1; 
the larger values give an increased rate of energy transfer and make the total 
coupling occur earlier. 
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These tests show that although the results can indeed be influenced by dif- 
ferent choices of the parameters, there is no serious qualitative change. We have 
also checked on the sensitivity of the code to the form chosen for the Eddington 
factor /^. (The background to our choice for this was discussed in Section II.) 
In our investigation, we have modified the form given in Eq. (5) by replacing the 
monotonic correction term X/{X + R) with a more sophisticated expression having 
a maximum whose position and the amplitude could be suitably tuned in different 
combinations. Paying attention to producing a smooth join between the optically- 
thin and optically-thick limits, we have found that reasonable variations in this 
joining function lead to only minor differences in the results, confirming previous 
experience [13, 14]. 

VI. Conclusion 

In this paper, we have presented a study of the hydrodynamics of the cos- 
mological quark-hadron transition in the presence of long-range energy and mo- 
mentum transfer by electromagnetically and weakly interacting particles. The 
relativistic radiative transfer problem for the generalized radiation fluid has been 
treated using the PSTF tensor formalism. A system of Lagrangian hydro dynamical 
equations has been presented which has then been solved numerically by means of 
a computer code which uses a standard Lagrangian flnite difference scheme for flow 
within the bulk of each phase together with a characteristic method for the vicinity 
of the phase interface across which relativistic junction conditions are solved. The 
results show the progressive coupling together of the strongly-interacting matter 
and the radiation fluid as the bubble expands. When the complete coupling oc- 
curs, there is no dramatic effect on the bubble which simply decreases its expansion 
velocity and the eventually approaches a similarity solution. 

Work is now in progress in order to extend this treatment to the flnal stages 
of the transition during which the evaporation of disconnected quark droplets oc- 
curs. At these stages, which are of great interest in connection with possible con- 
sequences of the transition, long-range energy and momentum transport certainly 
play an important role. 
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Appendix 

In this Appendix we outline the calculations by means of which the radiation 
characteristic equations (41)-(44) are rewritten in the new form (47)-(49). The 
aim is to make the ratio w\l jj^ the coefficient of a quantity which is small when 
/^ ^ and in order to do this we have first isolated the terms in (41) where this 
ratio appears: 
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This can be rearranged to give 
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Replacing (Al) by (A2) then leads to the expressions (47)-(49) which are the new 
characteristic equations for the radiation fluid. Note that the manipulation has 
produced two major effects: firstly it has introduced terms which are intrinsically 
small when /^ ^ (such as the difference between the sound speeds or /^ itself), 
and secondly it has brought together a set of terms which coincides with those 
on the left-hand-side of the standard fluid characteristic equation but are now 
evaluated along the characteristic directions of the radiation fluid. The sum of these 
terms becomes very small when /^ is small so that the two sets of characteristics 
nearly coincide. 
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Abstract 

Two previous papers in this series have presented a study of the growth of 
hadronic bubbles during the cosmological Quark-Hadron transition, treating 
the material within each phase as a single perfect fluid. Here, we extend the 
analysis to include the effects of long-range energy and momentum transfer 
by weakly and electromagnetically interacting particles. After a short review 
of the formalism adopted, we discuss the numerical strategies used in the 
computer code which has been constructed in order to solve this system 
of equations. Results for the growth of single hadronic bubbles are also 
presented. 
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I. Introduction 

At about 10 /xs after the big bang it is thought that strongly interacting 
matter underwent a transition from a plasma of free quarks and gluons to one in 
which the quarks were confined within hadrons (predominantly pions). Current 
lattice-gauge calculations [1] favour the view that this change may well have been 
a continuous one, but the possibility of it being a first order phase transition is 
by no means ruled out and considerable interest continues to be focussed on the 
astrophysical consequences which might have arisen if it was, indeed, first order 
(see, for example [2-4]). 

The work described here (which is part of an ongoing research programme [5- 
8]) is within the scenario of a first-order transition starting with the nucleation of 
hadronic bubbles in a slightly supercooled quark-gluon plasma. The bubbles then 
proceed to grow (with the quark-gluon plasma being progressively transformed into 
a hadronic one) until eventually they coalesce and give rise to disconnected quark 
regions which then shrink and probably disappear completely (possibly leaving 
behind a significant baryon-number inhomogeneity) [9] but might instead reach a 
stable configuration composed of strange quark matter [10]. 

In the earlier papers of the series, attention has been focussed on studying the 
hydrodynamics of the growth of single hadronic bubbles during the initial stages 
of the transition where it makes sense to consider the material in each phase as a 
perfect fluid composed only of the strongly interacting matter. While the transition 
involves only these particles in a direct way, an important role is also played by 
other particles present, which can interact with them through the electromagnetic 
and weak interactions: primarily photons, electrons, muons and their antiparticles 
(electromagnetic and weak interactions) and neutrinos and antineutrinos (weak 
interaction only). All of these have mean free paths long compared with that 
of the strongly interacting matter and can provide a mechanism for long-range 
transport of energy and momentum through the strongly interacting fluid. Since 
the leptons concerned are essentially massless, both they and the photons can be 
treated as components of a generalized "radiation fluid" and the problem is then 
one of relativistic radiative transfer. 

During bubble growth, the effect of this transport becomes significant when 
the radius Rs of the bubble surface becomes roughly comparable with the mean 
free path A of the particles concerned (~ 10'^ fermi for the electromagnetic in- 
teraction and ~ 1 cm for the weak interaction). When Rs <^ \, the bubble is 



essentially transparent to the radiation, which can then be ignored, while when 
-Rs ^ A the coupling is essentially complete on relevant length-scales so that the 
radiation and strongly interacting matter move together as a single fluid. Clearly, 
the process of coupling can, in principle, occur twice during the bubble growth, 
but because the behaviour is similar in each case, we will discuss here only the 
one occurring between the strongly interacting and electromagnetically interacting 
particles. Identical considerations apply also for the coupling with the neutrinos, 
the only difference, apart from the scale, being the different number of degrees of 
freedom into which the energy liberated by the transition is channeled. 

This paper is concerned with adding the effects of radiative transfer to our 
earlier scheme of calculation for following bubble growth [7], thus extending its 
range of validity. Subsequently, we aim to apply a similar scheme to the last stage 
of the transition (shrinking away of the disconnected quark regions) which is likely 
to give rise to the most interesting consequences from an astrophysical point of 
view. The formalism used for handling the radiative transfer has been the subject 
of an earlier paper [8] and only a brief outline of this will be presented again here. 

The plan of the paper is as follows. In Section II we give a short review 
of the equations and then Section III describes the outline of our computational 
scheme which uses a combined characteristic and Lagrangian flnite-difference ap- 
proach. Some particular problems arise in the computational implementation and 
we explain how these have been overcome by means of introducing new variables 
for the radiation fluid. Section IV concerns the setting up of initial conditions. 
Section V describes tests carried out and results obtained and Section VI is the 
conclusion. Throughout, we use a system of units in which c = % = k^ =1 and a 
space-like signature ( — ,+,+,+). Greek indices are taken to run from to 3 and 
Latin indices from 1 to 3; covariant derivatives are denoted with a semi-colon and 
partial derivatives with a comma. 

II. Relativistic Hydrodynamical Equations 

This section contains a short review of the equations used in the rest of the 
paper. A full discussion of these has already been presented in [8] to which the 
reader is referred for further details and references. 

The transfer of energy and momentum is considered as taking place between 
a standard fluid, composed of the strongly interacting particles and any other 



particles effectively moving together with them, and a radiation fluid consisting of 
those particles responsible for the long-range transfer of energy and momentum. 
We treat the transfer using the PSTF (Projected Symmetric Trace Free) tensor 
formalism of Thorne [11], [12] in which the relativistic radiative transfer equation 
is rewritten in terms of an infinite hierarchy of differential equations involving 
PSTF moments which is then truncated by introducing suitable closure relations. 
For spherical symmetry, the tensor formalism becomes effectively a scalar one. We 
made the truncation at second order and obtained 
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(1) 
W2 = asi, 

(2) 

W2 = Ie^O, (3) 



where wq is the energy density of the radiation, wi is the radiation flux and W2 is 
the anisotropy scalar of the radiation (all measured in the local rest frame of the 
standard fluid); sq and si are energy and momentum source functions and a and 
b are metric coefficients of the spherically symmetric line element 

ds^ = -a^dt^ + b^dfi^ + R^{de^ + sinH dip^), (4) 

where /x is a comoving radial coordinate and R is the associated Eulerean coordi- 
nate (the Schwarzschild circumference coordinate). 

The quantity f^ appearing in the closure relation (3) is a variable Eddington 
factor which can take values ranging from 0, for complete isotropy (which could 
be caused by the medium being extremely optically thick), to 2/3 for complete 
anisotropy (which might arise when the medium is very optically thin). The 
expression used for it has to be arrived at on the basis of physical considerations 
and for our case we have used 

^ 8u'/9 f A \ 

•^^ ~ (i + Auys) \\ + rJ' ^ ' 

where A is the effective mean-free-path of the radiation particles. In our picture, 
the bubble is first nucleated within a uniform and isotropic quark medium and the 
radiation field (which interacts with the quark medium on suitably large scales 
and is in thermal equilibrium with it) initially shares these properties [i.e. wq is 



constant everywhere and wi and W2 are zero). The radius of the newly-formed 
bubble is small compared with A and the medium is essentially transparent to the 
radiation on this scale. As the bubble starts to expand, the radiation quantities 
deviate from their values at the time of nucleation primarily as a result of the 
Doppler effect arising from the motion of the standard fluid rest frames with re- 
spect to that of the radiation field (which is, so far, remaining uniform in its own 
frame). These Doppler corrections can be calculated analytically and the results 
were presented in the Appendix of [7]. Solely on the basis of this consideration, 
one finds that f^ is given by (8/9)^^(1 + Au'^ /3)~^ . As the bubble grows to dimen- 
sions comparable with A, there is progressive coupling between the radiation and 
the standard fluid on the relevant length-scales and this interaction tends in the 
direction of making the radiation more isotropic as seen from the standard fluid. 
This effect is approximated by multiplying the Doppler term in (5) by a correction 
factor (in the large brackets) which has the effect of producing the right behaviour 
in the optically thin and optically thick limits and giving a physically plausible 
join between them. 

The appropriate value of A is not known accurately but, on the basis of 
elementary considerations, we have taken A sa 10'^ fm. Tests made in order to 
investigate the sensitivity of the numerical code to the values adopted for the 
various parameters will be discussed in Section V. 

For the source moments sq and si we use the expressions 

So = j{e-wo) + {so)^^ , (6) 

^1 = -^, (7) 

where {so)^^ is a term expressing the contribution due to scatterings and e is 
the energy density for radiation in thermal equilibrium with the standard fluid. 
Assuming that it roughly follows a black-body law, e can be written as 



30, r;, (8) 

with 5^ being the number of degrees of freedom of the radiation fluid and T^ the 
temperature of the standard fluid. 

Obtaining a suitable expression for (so)^^; is less straightforward. While 
detailed derivations have been made for simpler applications [13-14], the lack of 



precise knowledge about the interaction processes in the present case, has led us 
express {so)^^ by the simple absorption and emission term 

Msc = yCe-'u^o), (9) 

where a2 is an adjustable coefficient ranging between zero and one. Fortunately, 
the results of the numerical calculations turn out not to depend sensitively on the 
value chosen; a discussion of this will be given in Section V. 

Equations (l)-(3), describing the processes of transfer between the standard 
fluid and the radiation fluid, need to be solved together with the hydrodynamic 
equations 
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where p is the relative compression factor (which plays the same role as played 
by the rest-mass density in a non-relativistic fluid), T is the general relativistic 
analogue of the Lorentz factor, and w is the speciflc enthalpy (w = (e + p)/ p). 
The generalized mass function M can also be calculated using the alternative 
equation 

M^t = -47vR^R^t(p+-wo + -wi+W2]. (17) 

Equations of state are required for both phases of the strongly interacting matter. 
The hadronic medium is taken to consist of massless, pointlike pions for which 



while the quark phase is described by the M.I.T. Bag Model equation of state [15] 
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where B is the "Bag" constant. We take g = 37, §^ = 2 and these values need to 
be incremented by the relevant number of additional degrees of freedom 5^ when 
any of the radiative particles are completely coupled to the strongly interacting 
matter. 

For treating the phase interface, we again use a characteristic scheme (as in 
[7]). The characteristic form of (1), (2), (10) and (12) is 
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which are to be solved along the forward and backward characteristics of the 
standard fluid dp, = i: (a/b)cs dt (here Cg = {dp/deY'^ is the local sound speed in 
the standard fluid), and 
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which are to be solved along the forward and backward characteristics of the 
radiation fluid dfi = i: (a/fe)(l/3 + f^Y'^dt , and where for compactness we have 
defined 
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Supplementary equations calculating p, R and M are solved along the flow 
lines of the standard fluid (i.e. advective directions d/i = 0) and are 
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The configuration of characteristic curves adjacent to the interface is shown 
in Figure 1 for evolution of the system from time level i to a subsequent time level 
t + At. The dashed lines represent the forward and backward characteristics for 
the radiation fluid r, the full narrow lines are the equivalent characteristics for the 
standard fluids f, the vertical dotted line is a flow-line of the standard fluid in the 
quark phase and the heavy line is the worldline of the interface. 
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Figure 1. The configuration of characteristic curves near the phase interface drawn 
in the Lagrangian coordinate frame. 



Note that the difference between the characteristic directions results from 
the difference between the sound speeds in the radiation fluid ((1/3 + fsY ) ^^^ 
in the standard fluid (c^). If the latter were not relativistic, this difference would 



be large but in the present case Cg — > l/v3 and the difference between the sound 
speeds is frequently very small. This leads to some considerable complications in 
the numerical solution of the equations which we will discuss in Section IV. 

To complete the solution at the interface, it is also necessary to introduce 
junction conditions across it and those for the energy and momentum of the stan- 
dard fluid can be obtained using the Gauss-Codazzi formalism [16, 17]. Taking 
the surface tension a to be independent of temperature, these junction conditions 
are 
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where [A] = A'^ — A~ , {A} = A'^ + A~ , /x^ is the interface location, /i^ = 
dfig /dt, f = (a^ —b^fi^^Y'^ and the superscripts indicate quantities immediately 
ahead of and behind the interface [5] . 

Up to the time of the complete coupling, it is reasonable to neglect any 
interaction of the radiation fluid with the matter in the phase interface and so the 
energy and momentum junction conditions for the radiation are just continuity 
conditions: 
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Other supplementary junction conditions follow from continuity across the 
interface of the metric quantities R, dR/dt, and ds 
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At the time of bubble nucleation, conditions are essentially Newtonian so 
that 

[M]^=ATTRla. (34) 

One further equation is needed in order to complete the solution at the 
interface and this is an expression for the rate at which the quark matter passes 
across it. A suitable expression is obtained by setting the hydrodynamical flux Fj^ 
into the hadron region equal to the net thermal flux F^ into it 
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where ai is an accommodation coefficient (0 < ai < 1). 

This completes the set of equations. In the next Section we will present the 
details of the computational scheme and discuss the techniques and the strategies 
used. 

III. Numerical Computations of Bubble Growth 

In order to solve the equations of the previous section, we have constructed a 
computer code following the lines of the one developed previously for calculations 
of the initial stages of bubble growth (see [7] for a full description). As before, this 
employs a composite numerical strategy in which a standard Lagrangian finite- 
difference method is used to solve the hydrodynamical equations in the bulk of 
each phase, while the system of characteristic equations and junction conditions 
is solved for the grid zones immediately ahead of and behind the phase interface 
(which is tracked continuously through the grid). The radiation quantities wq and 
W2 are taken as "mid-zone" quantities while the radiation flux wi is taken as a 
zone boundary quantity and calculated at the half time level (as for the velocity 
u - see [7] for details). 

The equations presented in Section II are general in nature and can be applied 
to a variety of situations. Normally, there would be no problem in doing this but 
some particular difficulties have arisen when applying them to the present case of 
bubble growth at the cosmological Q-H phase transition. Here, direct use of the 
radiation equations in the form given above leads to rapidly-growing instabilities 
which destroy the solution. After a series of experiments it was found that the 



difficulty originates in the very small deviation of wq and wi away from their 
initial values during the early part of the bubble growth, and in the fact that 
the characteristic sound speed in the radiation fluid (1/3 + f^Y'^ becomes very 
close to the sound speed in the standard fluid c^ when the radiation is nearly 
isotropic in the rest frame of the standard fluid [i.e. when f^ — )• 0). These 
features lead to production of cancellation errors in the solution of Eqs. (1) and 
(2) and near divergences in the characteristic form of the equations (21) (the 
expression (c^ — 1/3 — ^ ) appears in the denominator of several terms). Note that 
the near equality of the sound speeds only arises when the standard fluid is, itself, 
relativistic (with Cg ~ l/v3)- Also, it is a peculiarity of the present situation that, 
initially, the radiation is nearly isotropic in the rest frame of the standard fluid 
not because the medium is optically thick on the scale of the bubble but, rather, 
because of the assumed isotropy of the universe. 

For overcoming the cancellation errors, we have introduced new radiation 
variables deflned as the difference between the energy density, the flux and the 
shear of the radiation fluid and some reference values (indicated below by the 
superscript *) with the aim of performing an analytic cancellation of large terms 
in the equations leaving behind smaller "difference" terms. It turns out to be 
convenient to take these reference values to be those which would be measured if 
the only effect were that resulting from the motion of the fluid relative to a uniform 
radiation fleld having an energy density equal to that at the time of nucleation 
of the bubble, [w(,)j^. (These are the pure Doppler values mentioned earlier and 
calculated in the Appendix of [7].) Using a tilde to denote the new variables, we 
have 



Wo = Wo - (wo)* = Wo - ( 1 + -u^ j {wo)j, , (36) 

4 
wi = wi - (wi)* = wi + -uT{wo)j^, (37) 

Q 
W2 =W2 - {W2y =W2 - -U^{wo)j,, (38) 



and equations (1) and (2) can then be rewritten as 
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Wo u \ M 
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SnaGRwi^ 



0, 



(40) 



where the partial derivatives with respect to /x have been replaced by the equivalent 
derivatives with respect to R [i.e. d/dR = [AttR'^ p/T)d/dfi). Equations (39) and 
(40) are the new radiation hydrodynamical equations for the bulk of each phase; 
once the "tilde" variables have been computed, the values of W(,,wi,W2 can be 
calculated from (36)-(38). Note that in (39), (40) the radiation variables which 
are multiplied by G are not transformed according to (36)-(38). This has been 
done to keep the expressions in a simpler form and because the contribution of 
these terms is small under the present circumstances. 

Using the new variables, the radiation characteristic equations become 



dwi±{- + f^ f'^dwQ + BU(Zm + ^Vdp + ^Tdt = 0, 



(41) 



where 



T.TT « J /'4 j,\ ~ 8, , [GM 



u^R 



±2(^ + /J^/^ 



2r 

u 
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3{\ + R) 
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(l+4ii2/3) 



(wo) 



1 + -u' 




(42) 
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The new form of the radiation junction conditions (28) and (29) is 
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(43) 



(44) 



«Vs ( 3 +/b j^o -(a^ +feVs)^i + i'^o) J, \abfi^ ll + -u^j ( 3 + /^ 



+ -uT{a'+b'fil) 



(45) 



a I 3 + /^ ) + ^ As r "^0 - 2ab^^wi 



+(wo), 



3 



a- I - + /^ + fe'/i 



2 .-.2 

s 



+ g^^^r/i, 



(46) 

Note that the characteristic equation no longer has terms with (c^ — 1/3 — f^) 
in the denominator but it does have terms containing the ratio wi/f^ and these 
still give rise to numerical instabilities. However, this can be countered by further 
rewriting the equations in a form in which wi/f^ only appears as the coefficient 
of expressions which are small when ^ — ^ 0. The central point in our strategy 
consists in isolating the group of terms which appears on the left-hand-side of the 
fluid characteristic equations (20) (and hence tends to zero when ^ — ^ and the 
fluid and radiation characteristics coincide). Details of the manipulation involved 
are given in the Appendix. This group of terms can then be conveniently handled 
using the differences between parameter values at the feet of the fluid and radiation 
characteristics. The revised form of the radiation characteristic equations (which 
is the one actually implemented in the code) is 



dwi±{- + f^ fl^dwQ + ^^]du - . ^^^ , dp + ^Tdt 
3 (e + p) 
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dt 



(47) 



13 




'R 



4 \ ^ 8, , fGM u 



± 



(1/3 + /J^/^ 



>'7V 



BT = (BU)<(G'<{47ri? 
2a 






Ml r 

+ ^ +— ^1 



+ 



4wi 



(48) 



P + ( 3 + /^ ) ^0 
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2u AtvGRwi 



-Tu{wo)j^ -wi)+ aBT, 



_ (49) 

where BT is the same as in (44). Note that the last two hnes of (47) are the terms 
on the left-hand-side of the standard fluid characteristic equations (20). 



IV. Initial conditions 

Making use of the experience gained in the earlier work [7], the setting of 
initial conditions has been quite straightforward. We started with a single super- 
cooled hadronic bubble nucleated in mechanical and thermal equilibrium with its 
surroundings at a temperature T^ slightly below the critical temperature for the 
transition T^ . The equilibrium is an unstable one and any perturbation (continu- 
ing expansion of the universe, for example) will cause it to start growing. However, 
this growth is extremely slow and, in practice, it is not easy to follow with our 
code as numerical noise rapidly dominates. Our strategy then was to introduce 
a small artificial perturbation, decreasing the fluid temperature inside the bubble 
by a small amount below its equilibrium value, and analytically tracing the effect 
of this on related quantities (including the velocity field of the standard fluid). 
This successfully produced a suitable data set for starting the time evolution in a 
smooth and consistent way (see [7]). 

Since, at this stage, the radius of the bubble is small compared with A, the 
radiation fluid is not significantly affected by the thermal perturbation and re- 
mains uniform and isotropic in its own frame. The initial conditions for W(,,wi,W2 
(measured in the rest frames of the standard fluid) are then those calculated from 
the Doppler formulae discussed earlier. 
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V. Tests and results 

As usual in numerical computations, the construction of the computer code 
was followed by a series of tests to eliminate errors and verify that the strategies 
used were satisfactory. One important test consisted in turning off the source 
functions and checking that the computed values of the radiation variables agreed 
with the analytical Doppler expressions. This revealed the problems discussed 
earlier. When these had been satisfactorily solved, the source functions were then 
turned on again and complete runs of the code were carried out. As the radius 
of the bubble increased (leading to increased coupling between the radiation and 
the standard fluid on relevant length scales) care was required as increasingly 
steep gradients of wq appeared in the vicinity of the interface prior to complete 
coupling. Since structure on a scale smaller than the grid spacing can obviously 
not be resolved, it is necessary to be ready to switch on complete coupling in the 
equations at the appropriate moment. Some experimentation was required in order 
to do this in the best way. When this had been done, further tests were carried 
out in order to examine the sensitivity of the results to changes in the physical 
parameters and assumptions. In the following, we will first present results for a 
set of "canonical" parameter values and then discuss the effect of varying some of 
these. 

In Figs. 2-5, we show results from a run with T^ = 150 MeV, T^ /T^ = 0.98, 
ctq = 1, ai = a2 = 1 and A = 10'^ fm. Here ao = cr/T^ is a parameter commonly 
used to measure the relative strength of the surface tension. The value which we 
are taking for this is larger than currently preferred ones but we give results for this 
case to allow direct comparison with those of [7]. Figs. 2 and 3 show the behaviour 
of the velocity and energy density of the standard fluid, at various times during 
the bubble growth, while Figs. 4 and 5 show the corresponding behaviour of the 
radiation energy density and flux. All of these figures should be viewed together 
and it is useful, also, to make comparison with Figs. 2 and 4 of [7] which are for 
the equivalent calculation without inclusion of the radiation particles. (Note that 
for convenience in drawing these figures, the values of the variables at the centre 
of the bubble have been plotted at /o5io-R(fm) = rather than at R = 0.) 

During the first part of the bubble expansion {^Rg ^ 10^ fm); the standard 
fluid variables behave in an identical way to that seen previously in the calcula- 
tion with no radiative transfer: the velocity of the interface progressively increases 
and a compression wave is pushed out into the surrounding quark medium. The 
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velocity profile in the quark phase is approximately solenoidal [u oc l/B?). The 
radiation variables at this stage have profiles which are almost exactly the Doppler 
ones produced by the motion of the fluid relative to an essentially uniform radia- 
tion field. In the previous calculations, where radiative transfer was not included, 
the standard fluid variables tended towards a similarity solution which was effec- 
tively attained for Rg ^ 10^ fm. In the present calculation, the coupling together 
of the radiation and the standard fluid (which can be clearly seen in Fig. 4) starts 
to be effective before the former similarity solution is fully reached, causing first a 
distortion of the velocity profile and subsequently a decrease in the peak velocity. 
(If a smaller value is used for (Tq, both bubble nucleation and the attainment of the 
similarity solution occur at smaller values of Rg-) As the coupling becomes more 
complete on the relevant length scales (when i^^ ~ 10'^ fm); the peak of the radi- 
ation flux profile becomes very narrow (the main part of the flux is concentrated 
exactly at the interface) and when it is no longer possible to resolve this on the 
grid we switch to total coupling. This involves setting to zero the radiation flux 
wi, the Eddington factor f^ and the source functions sq and si and augmenting 
the number of degrees of freedom for the standard fluid at the interface to include 
also those of the coupled radiation. The behaviour of the radiation at the inter- 
face is then included together with that of the standard fluid. Following the total 
coupling, the variables tend rapidly to a new similarity solution characterized by 
smaller velocities (the front now has to push a medium having larger inertia) and 
a smaller temperature jump across the interface. 




4 - 



2 - 



1 - 
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logio R (fm) logio R (fm) 

Figures 2-3. Velocity of the standard fluid u and energy density e. Different curves 
are for different times during the bubble growth; the dashed lines represent the 
values at the initial time. Here cco = 1 and A = 10 fm. 
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logio R (fm) 

Figures 4-5. Radiation fluid energy density Wq and energy flux w^. Here 0:2 = 1 
and A = 10 fm. 

Next, we turn to a discussion of the effect on the results of varying the values 

taken for some of the physical parameters. Bearing in mind the discussion already 

given in [7] , we will concentrate here just on the mean-free-path for the radiation 

particles A (whose value has an uncertainty of about an order of magnitude), and 

on the coefficient 02 appearing in the expression for the source moment sq (Eq. 

(9)). 
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Figures 6-7. Radiation energy density Wq and standard fluid velocity u just ahead 
of the phase interface and just behind it (higher and lower curves respectively). The 
curves are the result of calculations with different values of A; R^ is the bubble 
radius and 0:2 = 1. The curves are plotted so that the hadronic region is on the left 
of the discontinuity tracing the phase interface. 

The effect produced by varying A is illustrated in Figures 6 and 7, which show 
the behaviour during the bubble expansion of the values of wq and u measured just 
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ahead of the interface and just behind it. The different values of A used (which 
are represented with the different styles of line) are 5 X 10^ fm, 10'^ fm and 10^ fm, 
proceeding from left to right. As one would expect, use of a larger value of A has 
the effect of causing the coupling to occur when the bubble has reached a larger 
value of i?,. 
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Figures 8-9. Radiation energy density iuq and standard fluid velocity n just ahead of 
the phase interface and just behind it. The curves are the result of calculations with 
different values of the non-conservative scattering coefficient a.2- Here A = 10 fm. 



Figures 8 and 9 show the effect of varying the coefficient for the non- 
conservative scatterings (02). The different values of 02 used are 0, 0.5 and 1; 
the larger values give an increased rate of energy transfer and make the total 
coupling occur earlier. 

These tests show that although the results can indeed be influenced by dif- 
ferent choices of the parameters, there is no serious qualitative change. We have 
also checked on the sensitivity of the code to the form chosen for the Eddington 
factor f^ . (The background to our choice for this was discussed in Section II.) 
In our investigation, we have modified the form given in Eq. (5) by replacing the 
monotonic correction term A/(A + i?) with a more sophisticated expression having 
a maximum whose position and the amplitude could be suitably tuned in different 
combinations. Paying attention to producing a smooth join between the optically- 
thin and optically-thick limits, we have found that reasonable variations in this 
joining function lead to only minor differences in the results, confirming previous 
experience [13, 14]. 
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VI. Conclusion 

In this paper, we have presented a study of the hydrodynamics of the cos- 
mological quark-hadron transition in the presence of long-range energy and mo- 
mentum transfer by electromagnetically and weakly interacting particles. The 
relativistic radiative transfer problem for the generalized radiation fluid has been 
treated using the PSTF tensor formalism. A system of Lagrangian hydrodynamical 
equations has been presented which has then been solved numerically by means of 
a computer code which uses a standard Lagrangian finite difference scheme for flow 
within the bulk of each phase together with a characteristic method for the vicinity 
of the phase interface across which relativistic junction conditions are solved. The 
results show the progressive coupling together of the strongly-interacting matter 
and the radiation fluid as the bubble expands. When the complete coupling oc- 
curs, there is no dramatic effect on the bubble which simply decreases its expansion 
velocity and the eventually approaches a similarity solution. 

Work is now in progress in order to extend this treatment to the final stages 
of the transition during which the evaporation of disconnected quark droplets oc- 
curs. At these stages, which are of great interest in connection with possible con- 
sequences of the transition, long-range energy and momentum transport certainly 
play an important role. 
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Appendix 

In this Appendix we outline the calculations by means of which the radiation 
characteristic equations (41)-(44) are rewritten in the new form (47)-(49). The 
aim is to make the ratio wi/f^ the coefficient of a quantity which is small when 
_^ — )• and in order to do this we have first isolated the terms in (41) where this 
ratio appears: 
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This can be rearranged to give 
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Replacing (Al) by (A2) then leads to the expressions (47)-(49) which are the new 
characteristic equations for the radiation fiuid. Note that the manipulation has 
produced two major effects: firstly it has introduced terms which are intrinsically 
small when ^ — > (such as the difference between the sound speeds or f^ itself), 
and secondly it has brought together a set of terms which coincides with those 
on the left-hand-side of the standard fiuid characteristic equation but are now 
evaluated along the characteristic directions of the radiation fiuid. The sum of these 
terms becomes very small when f^ is small so that the two sets of characteristics 
nearly coincide. 
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